


THE 17" CHESAPEAKE SAILING YACHT SYMPOSIUM 
ANNAPOLIS, MARYLAND, MARCH 2005 


Toward numerical VPP with the full coupling of hydrodynamic and aerodynamic 


solvers for ACC yacht 


Erwan JACQUIN, Bassin d'essais des carénes, Val de Reuil, France, jacquin@bassin. fr 
Yann ROUX, Company K-Epsilon, Val de Reuil, France, yann@k-epsilon.fr 
Pierre-Emanuel GUILLERM, Bassin d'essais des carénes, Val de Reuil, France, guillerm@bassin. fr 


Bertrand ALESSANDRINI, Ecole Centrale, Nantes, France, bertrand.alessandrini 


ABSTRACT 


The classical approach of Velocity Prediction 
Program is to find the balance of hydrodynamic and 
aerodynamic sensors acting on the yacht to determine 
sailing conditions and associated performances. Usually, 
this approach is based on data given by _ towing tank, 
wind tunnel or numerical computations. We present in this 
paper the unsteady coupling between an hydrodynamic 
RANSE with free surface solver and a panel aerodynamic 
solver that allow to directly find the sailing condition of 
the yacht, and its performances, compute bay solving the 
six degrees of freedom motion equations of the ship in the 
hydrodynamic solver. The main advantage of this 
computational method is the decrease of numerical 
evaluation, compared to the classical interpolation 
approach to determine performances of a hull, and to 
directly rank several hulls in term of performances and not 
only in term of drag in fixed conditions. Further 
improvements will allow to simulate unsteady 
maneuvering of sailing yacht, and focus will for example 
on restart behavior after tacking, dynamic behavior in 
waves ... 


INTRODUCTION 


Based on results of hydrodynamic and aerodynamic 
studies, methods to predict performances of sailing yachts 
have been developed and integrated in software called 
Velocity Prediction Program (VPP). In 1936, Davidson 
first proposed a method to predict upwind performances of 
sailing yacht from towing tank data. A lot of improvements 
have been made in the development of those tools, 
especially by Murdey (1978), Herreshoff (1964), van 
Oossanen (1985), Oliver (1987) and other authors. When 
based on experimental or numerical results, the method 
used in the VPP always find the balance between 
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interpolated hydrodynamic and aerodynamic functions at 
full scale. Accurate results need to collect an important 
number of results from experiments, for example 200 
measurement points for one ACC hull during towing tank 
session. 


With the progress of numerical computation, Roux 
(2002) presented a new approach of VPP based on the 
coupling of aerodynamic and hydrodynamic computation. 
The hydrodynamic results were reported as interpolated 
surface function in the instationnary aerodynamic software 
which compute aerodynamic forces and find the balance of 
the boat. In the present paper we fixed the aerodynamic 
sensors and we compute the hydrodynamics sensors with a 
RANSE with free surface solver. At every time step we 
compute the hydrodynamic sensor and the aerodynamic 
ones to find, by the resolution of the six degrees of 
freedom motion equation the boat performances. 


The hydrodynamic and the aerodynamic solvers used for 
the simulations are first introduced, then we present several 
examples of unsteady hydrodynamic simulation of ACC 
yacht driven by an aerodynamic sensor. In order to show 
the ability of this method for direct hull performances 
evaluation, a final part presents the evaluation of two ACC 
hulls driven by the same aerodynamic sensor, in a 5 
degrees of freedom balance (yaw angle 1s fixed). 


HYDRODYNAMIC RANSE WITH FREE 
SURFACE SOLVER ICARE 


General presentation 


Numerical simulations of ship flows remain a 

great challenge in the hydrodynamics community: current 
effects, drift effects on a sailing boat, or manoeuvring 
effects for example, depend on the solution of this 
problem. The whole solution requires the computation of 
three dimensional unsteady turbulent boundary layers with 
flow separation connected to complex free surface effects 
and hull topology (appendages). 
In this paper we use Reynolds Averaged Navier-Stokes 
equations written under convective form in an unsteady 
curvilinear computational space fitted at each iteration to 
the hull and to the free surface. Fully non-linear free 
surface conditions are solved using an efficient fully 
coupled algorithm (Alessandrini & Al.) and turbulence 
effects are taken into account through classical k-@ 
turbulence model. 


Equations 
Governing equations, turbulence model 


The convective form of Reynolds Averaged 
Navier-Stokes Equations are written through partial 
transformation from Cartesian space (x/,x2,x3) to 


curvilinear space (€',é*,€°) fitted to the hull and the free 


surface at each time. Free surface elevation, the 3 Cartesian 
velocity components (ui), pressure (p) including 


gravitational effects (gx°) and turbulent kinetic energy 


2 
S pk ) are the dependant unknowns. 


Mean momentum transport equations are written in the 
moving referential : 
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where a; is the contravariant basis, g; the contravariant 


metric tensor, f; the control grid functions and u : the grid 


velocity which traduces the displacement of the mesh. 
Inertia forces due to non Galilean referential (gyration, 
accelerated translation) are taken into account in the gq; 
terms 


Mass conservation is expressed as the classical continuity 
equation : 


Finally, to close the equations set we used a classical 
k-—q@ turbulence model proposed by Wilcox [6], 


introducing a specific dissipation rate @ without low 
Reynolds formulation requirement. Transport equation of 
turbulent kinetic energy and dissipation rate are written as 
follows : 
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Free surface conditions 


Free surface boundary conditions are the 
kinematic condition, the two tangential dynamic conditions 
and the normal dynamic condition. Kinematic condition, 
coming from the continuity hypothesis, expresses that the 
fluid particles of free surface stay on it: 


h, + (bi (ub - ub hues} —¥ =0 


where b' is the bi-dimensional contravariant basis based 
on the discretization of free surface only. 


Dynamic conditions are given by the continuity of strains 
at the free surface. P pressure is assumed to be constant 
above free surface and normal dynamic condition is : 


where y is the surface tension coefficient (that is a 


physical way to smooth free surface near the hull) and r the 
free surface medium curvature radius. Tangential dynamic 
conditions are simply given by a linear combination of first 
order velocities derivatives : 


Peg ih az 
agg u_, =0 


Discretization 


General discretization method is based on second 
order (in space and time) implicit finite differences. 
Discrete unknowns are distributed on a_ structured 
curvilinear grid fitted to the hull and the free surface. 
Velocity Cartesian components, kinetic turbulent energy, 


and specific dissipation rate are located on the grid nodes. 
Pressure 1s located on the center of each volume and free 
surface elevation is located on the center of free surface 
interfaces. 


Convection terms are computed using an upwind 
second order scheme that needs a 13 nodes cell. Diffusion 
terms need 7 nodes for second order derivatives and 12 
nodes to express cross second order derivatives while 
pressure gradient requires 8 nodes for each component. 


Concerning the free surface calculation, it has 
been shown that the classical way using normal dynamic 
condition as Dirichlet condition on the pressure and 
uncoupled kinematic equation as transport equation to 
compute free surface elevation induces some problems 
connected to difficulties to exactly solve mass conservation 
under free surface. Efficient solution consists in using fully 
coupled algorithm (Alessandrini and Al.) that requires at 
each iteration the linear solution of mean momentum 
equations, continuity equation and all boundary conditions 
including free surface condition. Unfortunately the most 
efficient iterative algorithms (CGSTAB + ILU, Multigrid) 
are unable to invert this system because of the very bad 
matrix conditioning. The solution consists in modifying the 
system using free surface boundary conditions to express 
the flux through free surface. In this case, conditioning 
number decreases dramatically and fully coupled system 
can be inverted by iterative algorithms. 


Resulting linear system for velocity (U) and 


pseudo velocity (U ) components, pressure (P) and free 
surface elevation (#1) is written as follows 
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Example of validation on ACC yachts 


The first example of hull classification 1s 
presented for two ACC hulls. The two hulls are relatively 
close, and have been tested 1n tank. 





Figure 1 : Hull shape of ACC1 (left) and ACC2 (right). 


Figure 2 shows the computed and measured total 
drag at model scale on the ACC1 hull. The mesh size is 
composed of 300000 cells, and each speed is computed on 
approximately 10 hours on a single processor of a recent 
PC cluster. 
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Figure 2 : Comparison of computed and experimental 
results on the one ACC hull. 


The computed and measured results are very 
closed, but for ACC projects, the most important quality 
for a solver is to very accurately rank hulls in the same 
way than in tank. Figure 3 shows the total drag difference 
between the two hulls measured in tank and computed at 
model scale. The experimental results shows that the 2 
curves cross each other at the speed of 7.5 and 9.3 knots 
that are relatively well predicted by the computations. 
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Figure 3 : Total resistance difference in percents between 
two ACC computed with ICARE and measured in tank. 


On that example, the ability of the solver to rank 
hulls as in tank is shown. Even if the magnitude of the 
differences predicted by the solver are not exactly the 
same as those measured in tank, the ranking is always 
respected and all the cross of total drag curves are 
relatively well predicted. 





Figure 4 : Free surface elevation of ACCI (top) and 
ACC2 (bottom) 


AERODYNAMIC : PANEL METHOD SOLVER 
AVANTI 


General presentation 


The unsteady aerodynamic simulation is based on a 
description of the flow by means of lifting surfaces for the 
sails and particle method for their wakes. At each time 
step, particles are created along given lines (the trailing 
edge of the sail) in order to satisfy a condition formally 
derived from the Joukowsky condition. The velocity field 
of the particles 1s computed by the Biot-Savart relation and 
the deformation term is obtained through an integral 


relation by differentiating the Biot-Savart law. 


Equations 


Flow modeling is based on the vortex element 
method (VEM). This method is suitable for external flows 
for bounded vorticity support. It is the case for lifting 
surfaces, where the turbulent shear flow along the surface 
and the wake formed by the vortex shedding along the 
trailing edge are represented by dipole surface distribution 
and vortex sheets, respectively. During the last ten years, 
this method was successfully validated, particularly to 
capture dynamic wake effects in details. 


This method is basically made of two parts : a lifting 
body problem and a wake problem. These two problems 
are coupled using a kind of Kutta condition which has been 
derived from the kinematic and dynamic conditions along 
the separation lines. Usually, these lines are reduced to the 
trailing edges although more complicated situations have 
sometimes been considered. Except when writing this 
Kutta condition, the flow has been assumed to be inviscid. 
The lifting problem is solved using a boundary integral 
method : the surface of the body is represented using 
panels of rectangular shape which are used to satisfy the 
potential slip conditions. Specifically, a dipole strength 
was associated with each panel, and the strength of the 
dipole was adjusted by imposing that the normal velocity 
component at the surface of the body must vanish at 
control points (c.f. Fig.1; Eqs. 10~11). The wake has been 
modeled using the particle method itself (c.f. Rehbach, 
1977). According to this method, the vorticity distribution 
within the wake is described using particles carrying 
vortices. The particles motion is computed in a Lagragian 
framework. The vorticity on each particle has to satisfy the 
Helmholtz equation. 

According to the Helmholtz's decomposition theorem, the 
velocity field is: 


U=U,,+U,+U,+U,, 
(1) 


F is derived 


from a scalar potential ¢@ and U, , 1S derived from a 


where U, represents the inflow velocity, U 


potential vector yw, representing respectively the body's 


influence and the wake's influence, U., , is the velocity 


field induced for example by an another lifting body with 
its wake. 
The Euler equations in velocity-vorticity formulation for 


particle i are expressed in Lagrangian coordinates. X, is a 


point which can be the geometric center of the particle, the 
baricentric center for the vorticity contained 1n the particle 
or any other point representative of the particle location 


and 2. their circulation which is the amount of vorticity 


initially contained in the particle: 
Q. = (ly édo; O,rX,= (iy ONxdo. 
(2) 


where @ is the absolute vorticity. 


The evolution equation for X ie Q. are: 
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Let N; be the number of bound vortex particles equivalent 
to the dipoles (Hess, 1969), N,(t) the number of free vortex 
particles and N; the number of panels: 
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Figure 5 : Dipole-vortex equivalence 





The dipole distribution on the surface of the body, which 
is used to compute U g > has been set in order to satisfy the 


non-penetration condition. It is obtained by solving the full 
matrix linear system: 


[4 ][“]=[S] 
(6) 


where zis the unknown vector of the system, S the known 
vector of boundary conditions (-[0.. + U, + U,., | nN ] 


and A a square matrix called "influence matrix". A 
depends only on the body geometry. The a, term of this 
matrix represents the influence of the unknown ys; on the 
boundary term S;. It is obtained by computing the normal 
velocity induced by a unit dipole strength w at the control 
point of panel 7: 


(c.f. Eq. 11) 
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The singular behaviour of these equations when x, tend 


to _Y , leads to the introduction of a regularisation function. 
Thus, we have replaced the singular kernel (c.f. Eq. 10) by 
the convolution product of this kernel by f which 1s defined 
by: 


3 
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We set for the desingularised Biot-Savart law: 
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where O, is the distance with respect to the singularity, 


from which the regularisation is done. 
The amount of vorticity initially contained in the particle is 
obtained by application of the Bernoulli relation on the 


separation line (with U, _ =(o° 4 U) , the trailing edge 


velocity): 


Ou = — 
——+U, -sradu=0 
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So, during one step time of, we evaluate the location and 
the vorticity of the new vortex particle by the following 





equations: 
Q, =| 51, (us, (t+ At)— 4, (t)) ji + ar, Hate |j 
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(9) 


=> 


(x ee is the local frame of reference with nthe 


unit normal vector on panel i and i the tangent vector on 
the trailing edge. d/is the length of the side of panel 7 


which is a part of the trailing edge. 

For more generality, we consider two Cartesian frames of 
reference, the first one is fixed (inertial frame), and the 
other one moving relative to the first with the translation 


velocity V,, and the angular velocity OQ: 


=> 


U,(M)=V,(O)+®AOM 


The relation between the velocity U, (M) of a fluid 
particle in the absolute frame of reference and its relative 


al 


velocity U,(M) in the moving frame is: 


(M)-U,(M) 


The non-penetration condition takes into account the 
relative velocity. We have the same linear system that Eq. 
12 with the new known vector of boundary conditions: 


S(t) =—(U,, (1) +U, (1) +U,, (t)-U, (1) A 
The vorticity generation model which established the link 


between solid walls discretizing and the particle wake is 
the same, taking into account the body motion U ; (M ) 


el. _U, 
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The absolute vorticity (6 —rotU .) on each particle has to 


satisfy the Helmholtz equation written in the moving 
relative frame: 
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where dD — oe ae (U, grad | is the convective derivative 
Dt ot 


in the moving relative frame of reference. 


Example of validation on ACC sail plans 


Although the numerical method has already been 
tested for various aerodynamic problems, a set of 
validation tests has been performed for sails in order to 
gain confidence in the ability of the model to describe such 
flows. The results have been reported in Charvet (1992) 
and are summarised in Table 1. 


1.44 0.75} 1.10 
0.22} 0.066} 0.11 


1.37 0.80} 1.17 
0.13 0.05 | 0.10 


1.08} 1.37 
0.13 | 0.18 
1.08} 1.38 
0:095.|. 0.15 


Table 1 Comparison between computed and measured 
drag and lift coefficients on various jibs 


The agreement is quite good for the lift 
coefficient and acceptable for the drag one, considering 
assumptions made and uncertainties on the true shapes of 
the jibs. 








Figure 6 : Longitudinal and transverse forces on a ACC 
sail plane 
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Figure 7 : Wake around a ACC sail plane 


6 DEGREES OF FREEDOM MOTION OF THE SHIP 
General presentation 


The first step of the coupling of hydrodynamic 
and aerodynamic code is to solve the six degrees of 
freedom (6DOF) motion equations of the ship. This is 
done by solving the standard Euler’s law in the body fixed 
coordinate frame center on G. 


Coordinate frame system 


Navier-Stokes equations are solved is the fixed 
general axis center on 0 (Rj). This coordinate frame is 
Galilean, so acceleration terms of the fluid in the Navier- 
Stokes equations are not taken into account, and then 
reduce the complexity of the problem to solve. 

The Euler’s law are solved in the body fixed coordinate 
frame center on G (Rg). (R3) is an intermediate coordinate 
frame that is centered on G and in translation with respect 


to (Ro). 





Figure 8 : Fixed reference, horizontal body axis and body 
fixed axis use to solve 6DOF motion equation of the ship 





Figure 9 : Decomposition of the rotation of the ship 


The rotation matrix used to transform coordinates from the 
fixed coordinate frame to the body fixed coordinate frame 
is defined by : 


PRR me. Prior, PR, R, : 


And can be written : 


cly )e(0) 
Pr, or | cly)s()s()+ sy Je(g) 
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Motion equations 


The ship motion equations are written using the Euler’s 
law: 





dV, = & 2 
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dt 
d(ép) ee 
IR i +p -(TrOp)=Meyp 
Where 
m is the mass of the ship 
Ip is the inertial matrix of the ship 


Vo;pr 18 the velocity vector of the center of gravity 


is the angular velocity vector 
Fp is the total forces acting the ship (hydrodynamic, 


aerodynamic and gravitational forces) 
Mo;p_ 18 the total momentum acting the ship 


(hydrodynamic, aerodynamic and gravitational 
momentum) 


Those equations are discretized with explicit first 
or second order schemes. The use of explicit schemes 
means that at each time step, the hydrodynamic sensor is 
computed on the previous ship position. This leads to 
important instability that are solved, for example, using 
implicit or 4 order Runge-Kutta schemes. The use of 
theses schemes leads to important modifications of the 
hydrodynamic solver algorithm, and we prefer to use non- 
linear time step at each physical time step. That means that 
at each time step, the position of the ship is updated until it 
reaches a physical balance. 


Coupling with aerodynamic solver 


Classical computations with ICARE are carried 
with the moving speed, the heel and the leeway angle as 
input, with free pitch and trim. The first step of the 
coupling that is presented in this paper is to drive the 
sailing yacht with an external aerodynamic sensor pre- 
computed by the aerodynamic solver. 

As in the reality, the boat start without forward speed. 


Fx hyd, Fx Aero 


During few time steps, the speed is increased to the non 
dimensional speed of 1, with free pitch and trim. The boat 
is then driven with free surge. By solving the motion 
equations, the surge speed is adapted until the balance of 
all forces acting on the longitudinal direction 
hydrodynamic forces and aerodynamic driving force. 

The same procedure is used for free heel motion (including 
gravitational forces in the balance), and for sway motion in 
order to take into account the leeway on the hull. 
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Figure 10 : Time convergence of boat speed driven by a 
pure longitudinal aerodynamic force. 


EXAMPLE OF HYDRODYNAMIC COMPUTATION 
WITH EXTERNAL AERODYNAMIC DRIVING 
SENSOR 


In order to explain the interests of the sail driving 
computations, we present a succession of computations 
with several aerodynamic driving sensors acting on two 
ACC hulls. The first hull is a generic ACC hull, and the 
second is obtained by an increase of 15% the width 





Figure 11 : Front view of ACC and ACC+15% (obtained 


with a increase of 15% of the width) 


Longitudinal aerodynamic driving sensor 


The first example is a idealized downwind ACC, 
1.e. without any heel motion. In order to evaluate 
performances of one hull in a given wind condition, the 
classical approach is to compute the total drag for several 
boat speeds, and to find the intersection with the 
aerodynamic force (Figure 12). The total drag must then be 
computed for at least three speeds around the expected 
boat speed, and interpolated. 
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Hydrodynamic drag interpolated 


Aerodynamic force 
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Figure 12 : Sailing speed determination by VPP approach 


With the present approach, only one computation 
is required to find the boat speed. As explained in the 
previous section, the aerodynamic sensor is imposed as a 
driving force, and we let the sailing yacht find the balance 
of the sensors by changing the speed. This is done during 
the unsteady computation, solving the 6 DOF ship motion 
equations. (3DOF in this example : surge, pitch and trim) 


Figures 13 and 14 presents the boat speed 
convergence of the two ACC during unsteady 
computations. The main advantage is that only one 
computation is required to find the boat performances, 
against at least three for the hydrodynamic function 
interpolation. 
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Figure 13 : Unsteady convergence of boat speed. Solid 
line ACC, dashed line ACC+15%. Full range 
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Figure 14 : Unsteady convergence of boat speed. Solid 
line ACC, dashed line ACC+15%. Zoom 


The numerical results are summarized in the following 
table (non dimensional speed) : 


Table 2 : Converged boat speed (non dimensional) for 
ACC and ACC+15% with the same aerodynamical sensor 


This first example is for a single parameter computation 
(boat speed), pitch and trim being always free. In the next 
example, the heel is taken into account, with an 
aerodynamic side force imposed on the two hulls. 


Longitudinal and transverse aerodynamic driving 
sensor 


Let’s take academic example of a sailing yacht sailing 
upright without leeway : longitudinal speed, pitch, trim and 
heel are free. 


We want to find the performances of the hull (speed) in a 
specific wind condition, given by the driving and 
transverse aerodynamic force. If we want to use the 
classical method to find the boat speed, we need to built 
the hydrodynamic response surface of the hull, and 
computed and least three speeds, and three heel angle for 
each speed : at least 9 computations. (Figure 15) 
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Drag vs Boat speed for several heel angles 
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Figure 15 : Hydrodynamic response surface used for 


classical VPP approach (drag and lift vs boat speed for 
several heel angles). 


We can also make a single computation by 
imposing a driving and transverse aerodynamic force, and 
let the yacht find the balance of the hydrodynamic and 
aerodynamic sensors. 





Figure 16 : Free surface around the ACC and aerodynamic 
driving force represented by a vector. 


The solution is then the speed and heel angle of the yacht, 


with always free pitch and trim. 
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Figure 17 : Unsteady convergence of boat speed and heel 
angle. Solid line ACC, dashed line ACC+15%. 


The numerical results are summarized in the 
following table, which allow to compare the performances 
of the two ACC in a single computation (non dimensional 
speed) : 


| ACC 


Table 3 : Converged boat speed (non dimensional) and 
heel angle for ACC and ACC+15% with the same 
aerodynamic sensor. 





As we could expect, the wider hull as a lower heel 
angle, and lower speed. This due to the fact that the drag 
decrease with the heel, and that at a given heel angle, the 
drag of ACC+15%is more important than ACC one. 


Figures 18 and 19 shows the converged attitude of 
the two ACC hulls under the same aerodynamic sensor. 





Figure 18 : Free surface and heel angle of ACC sailing 
with the aerodynamic driving force. 





Figure 19 : Free surface and heel angle of ACC+15% 
sailing with the same aerodynamic driving force. (ACC 
boundary 1s plotted in red line) 


CONCLUSIONS 


Results presented in this paper show the ability to 
perform hydrodynamic unsteady computations with a 
RANSE with free surface solver, by imposing an 
aerodynamic sensor. This sensor is computed by a panel 
aerodynamic solver, using a particular method. Some 
results of such computations have been presented on two 
ACC hulls. The main interest of this method is to compute 
hull performances on the same aerodynamic sensor, on a 
single calculation. The results of the computation are then 
not the classical drag or lift of the hull, but boat speeds and 
motion (heel, leeway, pitch and trim). 


The next step of this work is to couple both 
solvers, in order to create a full numerical VPP, taking into 
account the full interaction of aerodynamic and 
hydrodynamic solvers. 
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